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Abstract Building circuits and studying their behavior in 
cells is a major goal of systems and synthetic biology. 
Synthetic biology enables the precise control of cellular 
states for systems studies, the discovery of novel parts, 
control strategies, and interactions for the design of ro- 
bust synthetic systems. To the best of our knowledge, 
there are no literature reports for the synthetic circuit 
construction for protozoan parasites. This paper describes 
the construction of genetic circuit for the targeted enzyme 
inositol phosphorylceramide synthase belonging to the 
protozoan parasite Leishmania. To explore the dynamic 
nature of the circuit designed, simulation was done fol- 
lowed by circuit validation by qualitative and quantitative 
approaches. The genetic circuit designed for inositol 
phosphorylceramide synthase (Biomodels Database — 
MODEL 1208030000) shows responsiveness, oscillatory 
and bistable behavior, together with intrinsic robustness. 
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Introduction 

In recent years, the disciplines of systems, synthetic, and 
chemical biology along with bioinformatics have gained 
importance as the embodiments of the future of biological 
sciences. Some of the greatest advances in biology and 
biotechnology are now arising at the intersection of the 
top-down systems approach and the bottom-up synthetic 
approach [1, 2]. Recent systems biology paradigms such 
as computational systems analysis, methods for quantifying 
time-dependent gene expression, and bioinformatics cata- 
loging of cellular parts enable synthetic biology. Interdepen- 
dence of systems and synthetic biology lies in the expansion 
and automation of the process of parts identification [3, 4]. 
A major goal of systems biology is to combine these parts to 
design biological circuits and study their behavior inside the 
cellular environment. The rapid capacity to sequence new 
genomes has sparked interest in equally rapid annotation 
capacities. As a result, databases of prokaryotic [5] and 
eukaryotic [6] motifs have been curated, enabling the anno- 
tation of different promoters [7], transcription factor binding 
sites, and terminators [8]. In addition to providing a repos- 
itory of natural biological parts, these tools allow synthetic 
biologists to design novel regulatory elements by combining 
motifs in many interesting ways [9]. 

Biological circuits are built by the assembly of biological 
parts [10]. A part can be described as the sequence of DNA 
that encodes certain information, associated with a specific 
function. Assembly of these parts (promoters, ORFs, ribo- 
some binding sites, etc.) give rise to a device, several such 
devices make up a system. With a myriad number of bio- 
molecules involved together with their dynamic behavior, 
there tends be a certain complexity involved within every 
system. Because of its parts-to-whole approach, a synthetic 
device is also comparable to an electronic circuit wherein 
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stimuli function as input and the final output is decided by 
the signals that modulate the behavior of the system [11]. 

Synthetic biology offers the potential to build constructs 
like engineered biomolecules, gene networks, and program- 
mable organisms to alter the mechanisms underlying the 
disease. The first genetic toggle switch was constructed in 
Escherichia coli by Gardner et al to understand the bistable 
and oscillatory behavior of the circuit [1]. Since then, syn- 
thetic biology has evolved right from building up of simple 
toggle switches to complex genetic circuits, replacing or 
rewiring them to achieve the desired cellular behavior. 

One of the many applications of synthetic biology is in the 
field of infectious diseases. This paper highlights the use of 
mathematical modeling with deterministic and stochastic 
approaches as a central component for the synthetic circuit 
design of inositol phosphorylceramide synthase (IPCS), an 
enzyme belonging to the sphingolipid metabolism of Leish- 
mania. IPC is present together with other sphingolipids and 
sterols in the membrane microdomains [12]. As IPCS is para- 
site specific and has no functional equivalent in mammals, it 
can serve as putative drug target fox Leishmaniasis [13]. 



Mathematical modeling serves as a critical link between 
the concept and realization of a biological circuit [14]. The 
constructed IPCS genetic circuit was modeled and simulated 
using diverse set of algorithms. This assumption was shown 
to greatly simplify the stochastic model thereby reducing the 
computational complexity. Once the circuit was constructed, 
validation was performed using the qualitative and quanti- 
tative approaches. The qualitative behavior of the circuit can 
be very informative in models with multiple steady states 
and showing switch-like or oscillatory behavior. It can be 
studied either over small parts of the parameter space by 
simply scanning over defined ranges of parameters and the 
initial conditions or by the global bifurcation analysis. On 
the other hand, qualitative analysis gives an idea as to which 
parameters offer the best success in achieving a desired 
behavior, whether a certain design can exhibit the desired 
function. Identifying the most crucial parameters depends 
not only on the mathematical analysis but also on the bio- 
logical feasibility. Characteristics such as strength promoter, 
transcript, and protein stability are quite variable. The changes 
in the characteristics of biological components can at best be 



Fig. 1 Genetic circuit for IPCS 
includes two genes (IPCS1 and 
2) forming the toggle switch, 
coupled to a repressilator 
comprising of three genes 
(Lad, TetR, and LambdaR) 
[Biomodels Database 
(MODEL1208030000)] 
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qualitative; however, it is also important to find parameter ranges 
that show robustness to variations. Model validity can be 
checked by comparison of the simulation runs with the quanti- 
tative experimental data such as the time course data or by the 
steady-state analysis [15, 16]. Based on the current approaches 
dealt, herein, we describe how different approaches of bioinfor- 
matics and systems biology could enable novel synthetic biolo- 
gy applications for the treatment of Leishmaniasis . 



Materials and methods 

Circuit design 

The genetic circuit was built using Tinker Cell, a Computer- 
Aided Design Software for synthetic biology (http://www. 
tinkercell.com). For the construction of the switch for IPCS, 
homologous sequences of Leishmania major were retrieved 
from GeneDB and Uniprot. Input sequences were in FASTA 
format or Genbank format. Parts were obtained from the 
Registry of Biological parts (http://partsregistry.org/cgi/ 
partsdb/Statistics.cgi), which is a collection of genetic parts 
that can be mixed and matched to build synthetic biology devices 



and systems. Grouping of the orthologous proteins was done by 
OrthoMCL. Selected parts were inserted for the construction of a 
model along with their DNA sequences. The circuit was built 
using two mutually repressible genes coding for IPCS that 
constitute the genetic toggle switch and a three repressible pro- 
moters that made up the repressilator. Transcription repression 
reactions were assigned for the two genes that are a part of the 
toggle switch as they mutually repress each other. Similarly, 
transcription repression reaction was assigned to each of the three 
genes that make up the repressilator [17]. Coupling between 
toggle switch and the repressilator was done by the insertion of 
the transcription repression reaction between the first gene of the 
toggle switch and first gene of the repressilator such that gene of 
the repressilator controls the gene of the toggle switch [18]. 

Parameters like protein degradation, promoter strength, 
transcription rate, dissociation constants, and the Hill's coef- 
ficient along with the reaction rates for the reaction were set. 
Parameter scan was done using the steady-state analysis and 
the genetic circuit was subjected to tau-leap stochastic simu- 
lation to study the behavior of the constructed circuit [19]. The 
logic gates were deciphered once the bistable and the oscilla- 
tory behavior of the genetic circuit was achieved. Truth tables 
obtained were converted into the circuit scheme via the 
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Fig. 2 Degradation graphs for all the proteins after the steady scan analysis, a Degradation graph for IPCS 1, b degradation graph for IPCS 2, c 
degradation graph for LacI, d degradation graph for Tet R 
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Case 1:Change in the Kd values of IPCS1(Kd =1 and 1.01) 
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Fig. 3 Graphs depicting the results for the stochastic simulation for the circuit by changing the Kd values of IPCS 1 



Karnaugh map (K-map), a method to simplify the Boolean 
algebra expressions (www.kamaugh.html). 

Model validation 

Bioconductor package based on the R programming lan- 
guage was used for the model validation which includes two 
packages — (a) Gene Regulatory Network Inference Using 
Time Series (GRENITS) and (b) Boolnet. The time series 
data required for the validation was obtained from Complex 
Pathway Simulator (COPASI; www.copasi.org) by using 
time course method for the stochastic simulation [20]. 



Qualitative network modeling was done by the Boolean 
method using BoolNet. BoolNet helps integrate methods 
for the synchronous, asynchronous and probabilistic Bool- 
ean networks. Apart from the reconstruction of networks 
from time series data, robustness analysis via perturbation 
and Markov chain simulations followed by the identification 
and visualization of attractors can be done [21]. Quantitative 
network modeling was done by the Bayesian approach 
using the GRENITS package which uses the time series 
data generated by COPASI. Network inference using ODE 
time series data was done which gives the probability of the 
genes included in the circuit, suggesting the regulators of the 



Case 2: Change in the Kd values of IPCS2 (Kd =1.01 and 1.02) 
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Fig. 4 Graphs depicting the results for the stochastic simulation for the circuit by changing the Kd values of IPCS2 
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Fig. 5 Logic gates representing the toggle switch (having two 
IMPLIES gates) and the repressilator (with three NOT gates) 



circuit along with the network uncertainty [22]. Pajek, a 
computer program for the analysis of large networks and 
visualization of the graphs was used to visualize the transi- 
tion states network [23]. 



Results and discussion 

Circuit design 

The genetic circuit for IPCS is composed of a toggle switch 
controlled by a repressilator. Genetic toggle switch is made 
up of two genes (IPCS1 and IPCS2) which code for repres- 
sor proteins and mutually repress each other, reflecting the 
bistable behavior of the toggle switch. The repressilator is 
made up of three genes namely lactose repressor (Lad), 
tetracycline repressor (TetR), and the lambda repressor 
(Phage lambda). The repressilator genes acts in cyclic man- 
ner with Lad repressing TetR, TetR repressing 

LamdaR and LamdaR in turn represses Lad, producing an 
oscillatory behavior. Coupling between genetic toggle switch 
and repressilator was developed between Lad and IPCS1 
where Lad regulates the IPCS1 and represses it (Fig. 1). 

Steady-state analysis 

For the steady-state analysis, default parameters were 
initially assigned to the genetic modules on the basis of their 



Table 1 Truth table for the ge- 
netic circuit, A and B are the 
input signals of the toggle switch 
and C represents the input signal 
of the repressilator 



Input Input Input Output 
(A) (B) (C) 



connectivity in the circuit and the Kd values for the genes were 
set to 1 (Electronic supplementary material (ESM) Table 1). 
The translation rates and transcription rates are set to 1 for all 
the genetic modules. Since every reaction is a transcriptional 
repression reaction, hill-hinze kinetics were applied to the 
genetic circuit and the value of hill's co-efficient was set to 
2; the kinetics have been included in the ESM. Steady-state 
scan analysis for the degradation rates of all the proteins was 
carried out. 

The degradation of IPC SI protein resulted in a rise in the 
expression of the opposite toggle gene IPCS2. At two time 
points, there was a peak fall in the level of IPCS2 after its 
degradation and here a switching behavior was seen (Fig. 2a). 
This desired graph was obtained at 0.18 degradation rate of 
IPCS1. Similarly on increasing the levels of IPCS2 at a 
particular interval of time, there was a decrease in the level 
of IPCS 1 as it was degraded, at these two time points, switch- 
ing behavior was observed (Fig. 2b). The desired behavior 
was obtained at 0.08 degradation rate of IPCS2. These two 
observations indicate the bistable behavior of the circuit 
wherein it was seen that IPCS1 is capable of repressing IPCS2 
and IPCS2 inrurn could repress IPCS1 and the flipping of the 
switch occurs. Initially, there was no coupling between Lac I 
and IPCS1 1, however at a time interval of 4 sec, there was an 
increase in the level of Lad along with a drop in the levels of 
IPCS1. This change in the regulation of IPCS 1 with respect to 
Lad shows the coupling between genetic toggle switch and 
the repressilator (Fig. 2c). In the degradation graph for TetR, 
there was a less oscillation between LambdaR and TetR; 
however, repression of LamdaR was observed (Fig. 2d). 

Simulation 

Simulation of the genetic circuit was performed for 100 time 
points using tau-leap stochastic simulation at differernt con- 
centrations. Based on the Kd values, the change in the 
protein level were studied, the results of fluctuations in the 
protein level with respect to change in dissociation constant 
(Kd) are dealt in two cases (cases 1 and 2). In case 1, it was 
observed that when the Kd value was 1, the level of IPCS 1 
was increased but when the value of Kd rose to 1.01, there 
was a switch in the level of IPCS2 as the dissociation 
between IPCS1 and regulatory binding site increased which 
lead the circuit to toggle increasing the levels of IPCS2. 
Similar fluctuations in protein level were also observed 
when the Kd values at regular intervals are altered (Fig. 3, 
ESM Fig. 11.2). Similarly in case 2, when Kd value was set 



Table 2 Karnaugh map 
of the truth table of the 
genetic circuit 
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Fig. 6 Digital circuit showing 
the coupling behavior of the 
toggle switch and the 
repressilator 




to 1.01, the level of IPCS2 was high but when the value of 
Kd rose to 1.02, a switch in the level of IPCS 1 was observed 
as the dissociation between IPCS2 and repressor increased. 
This caused the circuit to flip and level of IPCS 1 increased. 
Similar perturbations in the levels of IPCS 1 and IPCS2 were 
also observed when there was a change in the Kd at regular 
intervals (Fig. 4, ESM Fig. 11.3). These observations indi- 
cate that there is flipping of the toggle switch indicating that 
the constructed toggle switch comprising of IPCS 1 and 2 
exhibits bistable behavior and that there was a direct cou- 
pling between IPCS1 of the toggle switch and Lad of the 
repressilator. 

Digital circuit representation of the synthetic genetic circuit 



Digital circuits gives valuable insight into the logic gates 
implemented in the circuit design to understand the func- 
tioning of the circuit. In comparison to an electronic switch, 
the two mutually repressible genes constituting the toggle 
switch are comparable to two IMPLIES gates such that 
output of one IMPLIES gate is the input signal for the 
opposite IMPLIES gate. On the other hand, the repressilator 
is constructed by integrating an odd number of NOT gates in 
a circular fashion such that the output of the one is the input 
of another and the output of the last gate serves as an input 
of the first one (Fig. 5). Odd number of NOT gates are 
included in the repressilator such that the last output signal 
oscillates between high and low alternatively. 

The logic behind the IPCS genetic circuit can be repre- 
sented by a truth table (Table 1). Karnaugh map for the 
circuit was generated which is a simplified representation 
of the truth table that allows the conversion of the truth table 
into a Boolean formula (Table 2). This method derives two 
different descriptions of every digital circuit that are called 
as product of sums (POS) and sum of products (SOP). 
Digital circuit shows the logic gate for the toggle switch 
and repressilator with the SOP=(|,4 xQ+(|ix \B) and POS= | 
A*(\B+C). The circuit scheme for the circuit is represented 
in Fig. 6. 



Genetic circuit model validation 

Model validation was done by using qualitative and quanti- 
tative approaches. ODE model for the genetic circuit was 
generated giving an insight into regulatory mechanism for 
the circuit. The equations for the construction of ODE 
model are included in the ESM. 

Network inference 

Gene regulatory network was constructed for the genetic 
circuit using the Bioconductor package. The model was 
simulated using COPASI to obtain the time series data for 
qualitative and quantitative network modeling. Linear net- 
work construction was done as the circuit model was found 
to be stable. Qualitative network modeling was done using 
the GRENITS package which gave the probability of each 
regulator in the regulatory network circuit. Monte Carlo- 
Markov chain simulation was performed using the default 
parameters to obtain the posterior probability (Fig. 8). Dur- 
ing the simulation, two Markov chains are generated and 
based on their convergence the link probability of the net- 
work is deduced. A probability matrix for each gene in the 

Table 3 Posterior probabilities for each network connection 



From 



To 



Probability 



IPCS_1 

IPCS_1 

IPCS_1 

IPCS_1 

IPCS_1 

lad 

lad 

lad 

LamdaR 

LamdaR 



IPCS_1 

IPCS_2 

Tetr 

lad 

LamdaR 

IPCS_1 

Tetr 

lamdaR 

IPCS_1 

IPCS 2 
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Table 4 Posterior probability 
for number of regulator for each 
gene, 1 represents regulation 
between the regulators and 
indicates no regulation 
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regulatory network, an analysis and a convergence plot 
was generated for the network circuit. Network infer- 
ence was made for 10 and 100 time points, the network 
inference was observed to understand the regulatory 
mechanism in the circuit wherein a probability of 1 
shows the regulation between the respective regulators, 
while a probability of indicates no regulation. A 
probability of 1 was obtained between IPCS1 and 2 



and they appear to be regulating each other; in addition, 
a probability of 1 was also seen between Lad and IPCS 
1 reflecting the coupling mechanism (Tables 3 and 4). 

Analysis plot of the genetic circuit shows the link prob- 
ability which provides insight into the switching behavior, 
i.e., between the Lad belonging to the repressilator and 
IPCS1 of genetic toggle switch. A probability of 1 is repre- 
sented by blue color in the plot. The analysis plot shows the 
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Fig. 7 a Heat map plot of Network Link probabilities, b network uncertainty plot, c circadian clock network showing the switching behavior of 
IPCS1 andIPCS2 
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probability of 1 between the two genes IPCS1 and IPCS2. It 
was also observed that the LacI and IPCS1 genes appear to 
be linked as the probability between them was also 1 
(Fig. 7a). Network uncertainty is given by the marginal 
network uncertainty plot giving and idea about the top 
regulators in the circuit. From the network uncertainty plot, 
IPCS1, IPCS 2 together with LacI were the key regulators in 
the network model (Fig. 7b, c). 

Convergence plots The posterior means of the variables 
(gamma, B, lambda, and mu) each are compared in Fig. 8. 

Inferred network The inferred network displays the regula- 
tion between IPCS1 and IPCS2 which reflects the switching 



behavior of the toggle switch. Network inference made at 
10s time point shows no regulation between LacI and IPCS1 
and there was no coupling seen between the toggle switch 
and repressilator. It was observed that IPCS2 repressed the 
gene expression of IPCS 1 while LacI and TetR regulated the 
expression of LambdaR. IPCS1 remained in an ON state 
while IPCS2 was in an OFF state (Fig. 9a). However, at the 
second time point, i.e., network inference made at 100s, 
coupling between the toggle switch and the repressilator 
was seen which resulted into a repression of IPCS 1 by LacI, 
with IPCS1 in OFF state and IPCS2 in ON state (Fig. 9b). 
Thus, the probability derived justifies the associated param- 
eters that accounts for the designability of the IPCS genetic 
circuit. 



Gammas 




Fig. 8 a Indicator variable of Gibbs variable selection, b coefficient of linear regression, c precision of each regression, d intercept of each 
regression 
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Inferred network of the genetic circuit of IPCS at 10 sec 
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Fig. 9 Inferred network of the genetic circuit of IPCS showing the coupling between the genetic toggle switch and the repressilator 
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Table 5 Transition of states in 
the Boolean network for IPCS 



Boolean network with five genes involved 
Gene 1, gene 2, gene 3, gene 4, gene 5 
Transition functions 

Gene l=</(gene 3, gene 2, gene 4, gene 1, gene 5){ 1 1 101 101 101 1 1 101 1 1 1 1001 1 1 1 1 10000}> 
Gene 2=</(gene 3, gene 4, gene 5, gene 2, gene 1){00100001000010110100101000101011}> 
Gene 3=</(gene 5, gene 3, gene 1, gene 2, gene 4){10101110001010011000001111000110}> 
Gene 4=</(gene 4, gene 2, gene 3, gene 5, gene 1){1 1 10101001 1 101 1 1001 11000100001 1 1}> 
Gene 5=</(gene 3, gene 1, gene 5, gene 2, gene 4){0001001 1000101 1 1 1 1 10000001 1 1 101 1}> 



Qualitative network modeling 

To test the robustness of the circuit and identification of the 
number of attractors that circuit exhibits, qualitative Bool- 
ean network modeling for the genetic circuit was performed. 
After obtaining the steady state of the network, network was 



perturbed to check the robustness of the circuit. Boolean 
network for the circuit was constructed using time series 
data generated by COPASI. For the construction of Boolean 
network, binarization of the time series data was performed. 
The table containing the binarized time series data is includ- 
ed in the ESM. 
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Fig. 10 a Attractor plot for the initial attractor of the network, b Attractor plot for two steady states, c Sequences of path to attractor 
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Fig. 11 Robustness plot of the genetic circuit 
Boolean network 

Boolean network was constructed from the time series data 
by using the best- fit algorithm. The number of states present 
in the Boolean network generated for IPCS was 32. These 
32 states are represented by 1 (active state) and (inactive 
state). State transition can be observed, after which the 
network circuit attains steady states (attractors). Genetic 
circuit had two attractors that indicate the bistability of the 
circuit constructed. The two steady states obtained represent 
the ON and OFF state for each gene in the circuit, speci- 
fies the OFF state and 1 specifies the ON state. Attractor 2 is 
a simple attractor consisting of two state(s) and has a basin 
of 30 state(s). Genes are encoded in the following order: 
gene 1, gene 2, gene 3, gene 4, and gene 5. Transition of 



different states to obtain the steady state is represented by the 
transition table of the states which also includes the probabil- 
ity for each transition of the state (Table 5, ESM Table 1). The 
last state, represents the OFF state for IPCS1 and ON state 
for IPCS2 and ON and OFF state for the other genes, respec- 
tively. This reveals the switching of genetic toggle switch of 
circuit where IPCS1 was repressed by IPCS2. 

Attractors of the circuit are represented by the attractor 
plot which includes ON and OFF state for each gene in the 
network. In initial state, IPCS1 is ON, i.e., active state and 
IPCS2 is OFF, i.e., inactive state (Fig. 10a). Lad was in an 
OFF state, there was no coupling between Lad and IPCS1 
resulting in repression of IPCS 2 by IPCS1. After the tran- 
sition of the states, two steady states that were obtained are 
represented by an attractor plot which displays the two 
steady states justifying the dynamic behavior of the circuit 
to act as bistable genetic switch. The percentage to attain the 
steady states of the circuit was 93.75 % which reflects the 
robustness of the attractors (Fig. 10b). Attractor plot repre- 
sents the OFF state of IPCS 1 and ON state of IPCS2 which 
depicts the switching behavior of the circuit as IPCS2 
represses IPCS1 which was not in observed in the initial 
state. This shows that the circuit toggles periodically be- 
tween the two states and the two genes appear to repress 
each other mutually. This is in line with the hypothesis laid 
with the constructed genetic circuit for IPCS. Visualization 
of a sequence of states can be represented by the path of 
attractor 's plot where the columns in the table represent 
consecutive states of the time series. The last state is the 




Fig. 12 State graph of the transition of states (visualized using Kamada-Kawai layout) 
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steady-state attractor of the network. This gives the inactive 
and active state of each gene at different time points 
(Fig. 10c). 

Robustness of the genetic circuit 

Perturbations in the genetic circuit are introduced to test the 
robustness of the networks to noise and mismeasurements. 
The robustness plot for the genetic circuit and the percent- 
age obtained is greater than 50 % confirming the robustness 
of the constructed IPCS genetic circuit (Fig. 11). The net- 
work for the genetic circuit was exported to Pajek for the 
visualization of the transition states in the network (Fig. 12). 
Thus, the robustness plot of the circuit and the attractors 
obtained for the genetic circuit justifies the design of the 
IPCS circuit for Leishmania. 



Conclusion 

Genetic circuit designed for IPCS was subjected to simula- 
tion and validated by both qualitative and quantitative 
approaches. The circuit exhibited both bistable and oscilla- 
tory behavior along with robustness to perturbations. The 
methodology incorporated in this paper specifies those 
regions where the probability of observing the desired be- 
havior is appreciable. This may allow a more comparative 
assessment of different design proposals for other protozoan 
parasites, especially when dynamics are expected or desired 
to exhibit elements of stochasticity. Collectively, the advan- 
ces in the field of synthetic biology with the help of bio- 
informatics will help mitigate the great complexity and 
uncertainty currently impeding the study and design of 
synthetic circuits embedded in cellular systems. Such prog- 
ress holds the promise of streamlining and intensifying the 
whole drug discovery process chain and aid delivery circuits 
as next generation therapeutics. 
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